Initialisation

## [1] "number of cores available = 1"
#Phi[1] ; eta = valeur de fin  Phi[2] = valeur du noeud  Phi[3] = echelle
m <- function(t, eta, phi) (phi[,1] + eta)/(1+exp((phi[,2]-t)/phi[,3]))
#=======================================#
param <- list(sigma2 = 0.05,
              rho2 = 0.1,
              mu = c(5,4,0.5),
              omega2 = c(0.5,0.1,0.01))
#=======================================#
t <- seq(2,6, length.out = 10) #value of times

# --- longitudinal data --- #
data <- NLME_data(G = 40, ng = 12, time = t, fct = m, param = param)

getDim(data)
##    G   ng    n    N   F. 
##   40   12 4800  480    3
Y <- data$obs

SAEM avec simulation par MCMC

\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]

Phi <- fct_vector(function(sigma2, rho2, mu, omega2) 
  c(-n/(2*sigma2), -N/(2*rho2), -G/(2*omega2), G*mu/omega2), dim = c(1,1,F.,F.) )$eval

S <- fct_vector(function(eta, phi) mean((Y - get_obs(data, eta = eta, phi = phi) )^2 ),
                function(eta, ...) mean(eta^2),
                function(phi, ...) apply(phi^2, 2, mean),
                function(phi, ...) apply(phi  , 2, mean), dim = c(1,1,F.,F.) )

Metropolis Hastings

loglik.phi <- function(x, eta, Phi)
{
  id <- c(1,3,4) #indice de S_1 et S_{3,.} puis S_{4,.}
  sum( Phi%a%id * S$eval(eta, x, i = id) )
}
loglik.eta <- function(x, phi, Phi) 
{ 
  id <- c(1,2)
  sum( Phi%a%id * S$eval(x, phi, i = id) )
}

SAEM

Initialisation

# ---  Initialisation des paramètres --- #
para <- param %>% lapply(function(x) x* runif(1, 1,1.2))
para$rho2 = 0.1 ; para$mu <- c(4,5,2) ; para$omega2 <- rep(.1,3)

# --- Initialisation des chaines MC : Z_0 ---
Z <- list(eta = list( rnorm(G*ng, 0, para$rho2)  %>% matrix(ncol = 1) ),
          phi = list( matrix(rnorm(F.*G, para$mu, para$omega2), nrow = F.) %>% t ) )

Etape simulation et maximisation du SEAM

sim <- function(niter, h, Phih, eta, phi)
{
  M <- length(phi)

  eta <- 1:M %>% lapply( function(i)
    MH_High_Dim_para_future(niter, eta[[i]], sd.eta(h) , loglik.eta, phi[[i]], Phih, cores = 1))
  
  phi <- 1:M %>% lapply( function(i)
    MH_High_Dim_para_future(niter, phi[[i]], sd.phi(h), loglik.phi, eta[[i]], Phih, cores = 1 ))
  
  list(eta = eta, phi = phi)
}

maxi <- function(S)
{
  list(sigma2 = S%a%1,
       rho2 =   S%a%2,
       mu =     S%a%4,
       omega2 = S%a%3 - (S%a%4)^2 )
}
sigma2 rho2 mu1 mu2 mu3 omega21 omega22 omega23
Oracle 0.0494145 0.0989822 5.135294 3.995671 0.4979819 0.3929559 0.1507572 0.0111813
Initialisation 0.0511214 0.1000000 4.000000 5.000000 2.0000000 0.1000000 0.1000000 0.1000000 
niter <- 25
correction.phase <- 50
MH.iter <- function(k) ifelse(k<correction.phase, 20, 10)

sd.eta <- function(k) if(k<correction.phase) .03 else .04
sd.phi <- function(k) if(k<correction.phase) c(.028, .036, .021) else .01

res <- SAEM(niter, MH.iter, para, Phi, S$eval, Z, sim, maxi, 
            burnin = 150, coef.burnin = 3/4, 
            eps = 1e-3, verbatim = 2)

saveRDS(res, 'saem.rds')
gg <- plot(res, true.value = param)
## [1] "SAEM execution time = 02min 47sec"
Result of the SAEM-MCMC
sigma2 rho2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3
Real value 0.0500 0.1000 5.0000 4.0000 0.5000 0.5000 0.1000 0.0100
Estimated value 3.3205 0.1738 3.8711 5.0436 2.0748 0.2983 0.4205 0.1780
Rrmse 65.4090 0.7377 0.2258 0.2609 3.1495 0.4034 3.2051 16.7954